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Abstract 

Using equilibrium and non-equilibrium molecular dynamic (MD) simulations, we found that en- 
gineering the strain on the graphene planes forming a channel can drastically change the interfacial 
friction of water transport through it. There is a sixfold change of interfacial friction stress when 
the strain changes from —10% to 10%. Stretching the graphene walls increases the interfacial 
shear stress, while compressing the graphene walls reduces it. Detailed analysis of the molecular 
structure reveals the essential roles of the interfacial potential energy barrier and the structural 
commensurateness between the solid walls and the first water layer. Our results suggest that the 
strain engineering is an effective way of controlling the water transport inside nano-channels. The 
resulting quantitative relations between shear stress and slip velocity and the understanding of the 
molecular mechanisms will be invaluable in designing graphene nano-channel devices. 
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I. INTRODUCTION 

Graphene, due to its extremely large specific surface area superior electronic and 
mechanical properties, is a promising material in the novel nano-fiuidic devices applica- 
tions, such as, water desalination, nano-filtration, photo-catalysis, super-capacitor etc. 
0]. Graphene layers usually self-assemble in paper-like structures with interlayer distance 
on the nanometer scale 0, 3]. Fluid transport in such flat nano-channels is a key concern in 
improving the performance of the nano-fiuidic devices. Water confined inside nano-channels 
exhibits superior transport properties owing to its having a different molecular structure 
than the bulk water 8|, |9| . Experiments have shown that the flow rate of water inside car- 
bon nanotubes (CNT) is several orders of magnitude higher than predicted by conventional 



hydro dynamics 



10Nl2|. Because of their structural similarity, the extremely fast water 



transport found inside CNT is also expected in graphene nano-channels. Deep understand- 
ing of the molecular transport and the means of controlling the water transport inside the 
graphene fiat nano-channels are highly desired in practice. 

The structural properties of the single atomic layer of graphene, to a large extent, are 
determined by its interactions with the environment. For example, an epitaxial strain of 



±1% builds up in the graphene when it is grown on different substrates 13l-ll6|: elec 



tromechanical strain on the order of ±1 — 2% results from charge injection in the graphene 



171]; and applying mechanical force (e.g., nano-indentation) on the graphene can lead to a 



strain of about 10% [18]. Meanwhile, due to the size confinement in the CNT and graphene 
nano-channels, the first water layer next to the solid walls tends to be highly ordered 19| . It 
is reasonable to expect that the slip flow of water over the graphene layers will be analogous 
to the friction between two ordered crystal planes, meaning the interfacial commensurate- 
ness will play an essential role. Strain engineering of graphene layers to alter the interfacial 
molecular structures could thus provide a method of controlling the flow in nano-fluidic 
devices. 

Here, we use both equilibrium and non-equilibrium MD simulations to study the water 
flowing in uniformly biaxial strained graphene nano-channels. The relations between the 
interfacial friction shear stress r and the slip velocity Vg are obtained. Our results show that 
the interfacial friction coefficient (e.g., ratio of r over Vs) increases almost 6 fold when the 
strains applied to graphenes varies from —10% to 10%. To gain a deep understanding of the 
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physics behind this remarkable change, the molecular mechanisms, such as the interfacial 
potential barrier, density and structure factor of the first water layer, are quantitatively 
investigated. 



II. SIMULATION DETAILS 

Both equilibrium and non-equilibrium MD simulations were performed. Our molecular 
system is illustrated in Fig. dJ^a): water flows between two graphene sheets with an interlayer 
distance about 2 nm. Uniform in-plane biaxial strain is applied ranging from —10% up 
to 10% to the graphene walls. Accordingly the carbon-carbon bond length acc changes 
from O.Qacco to l.lacco; where acco denotes the carbon-carbon bond length of strain-free 
graphene (e.g., acCo = 0.142nm). To ensure the internal pressure of water at latm, we fix 
one of the graphene walls and use the other one as a piston to impose the pressure at a 
given temperature (300K). The position of graphene walls is then fixed in the following MD 
simulations. Periodic boundary conditions are applied along the two directions in the plane. 
Along the flow direction (x direction in Fig. [U^a)), the simulation box length was chosen to 
be 6.0 nm. However, we also performed simulations with a longer length 20 nm for the zero 
strain case and found the differences in the calculated friction stress were negligible. There 
are about 1200 to 2000 water molecules in our systems. 



All MD simulations were performed with the LAMMPS code 20|. A time step of 1.0 



fs was used and the total simulations time was about a few nanoseconds. We used the 



CHARMM force field and the SPC/E model |2lj for water with the SHAKE algorithm |22|. 



The water-carbon interactions were described by a Lennard- Jones potential between oxygen 
and carbon atoms with parameters e = 4.063 meV and a = 0.3190 nm, yielding a contact 



angle of 95° between the water and graphene 23| . The van der Waals forces were truncated at 



umbic interactions computed using the particle-particle particle- 



1.2 nm with long-range Co^ 
mesh (PPPM) algorithm 2^. Water molecules were kept at a constant temperature of 
300 K using the Berendsen thermostat, with the temperature calculated after removing the 
center-of-mass velocity. We also tried the Nose-Hoover thermostat applied to the degrees of 
freedom perpendicular to the flow direction in the zero-strain case. No significant differences 
between these two thermostats were found in our study. 

In our non-equilibrium MD simulations (NEMD), the Poiseuille flow was driven by ap- 
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FIG. 1: (color online) (a) Water flow inside graphene flat nano-channel with height about 2 nm. 
(b) The Plug- like velocity profile of the Poiseuille flow along the flow direction (i.e., x direction) 
observed in our MD simulations. We applied three different accelerations to water molecules to 
obtain the different steady flow velocities as shown by the different symbols. 

plying a constant gravity-like acceleration to all the oxygen and hydrogen atoms. Different 
external accelerations were applied to achieve a set of steady flows of different velocity (once 
the external forces were balanced by the friction). It usually took a few hundred picoseconds 
to reach a steady flow and the simulations were then continued for five more nanoseconds 
to collect data. By applying different accelerations from to 0.004 nm/ps^, we obtained 
the steady-flow velocities from to 200 m/s. Velocity profiles of the Poiseuille flows in our 
nano-channels are plug-like as shown in Fig. [D^b), so we can set the slip velocity equal to 



the average velocity. We noted that the similar simpliflcation was adopted in Ref jsj. The 
friction shear stress was calculated from the external force as r = Nma/2A, where N is the 
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number of water molecules, m is the mass of one water molecule, a is the acceleration, and 
A is the area of one graphene sheet. For comparison, we also calculated the shear stress r 
using the forces in the flow direction acting on carbon atoms (due to pair-wise interactions 
between the carbon atoms and the water molecules) and found good agreement, e.g., the 
difference is less than 5%. 



III. RESULTS AND DISCUSSIONS 



A. Strain dependent friction coefficient from MD simulations 

Figure EJ^a) shows the interfacial friction shear stress r as a function of flow rate Vg in 
the graphene channels placed under different strain. They are fitted well using an inverse 
hyperbolic sine (IHS) relationship t/tq = asinh(fs/fo). The values of tq and vq derived from 
the fits are summarized in table [T] for different strain values. The IHS relation's derivation 
arises from the transition state theory model of Yang [23] , in which the slip flow is described 
as a collective thermal diffusion of fluid atoms over a periodic solid wall. We have previously 
justified the use of this relation by using large-scale NEMD simulations for water flowing 



inside CNTs 
of strain. 



91]. Figure [2]^a) shows that the friction shear stress increases with increasing 



TABLE L Fitted parameters tq and vq in the inverse hyperbolic-sine relationship, t/tq = 
asmh(vs/vQ), which describes the water slip flow in the strained graphene nano-channels. The 
friction coefficient A and slip length Is are also calculated (see text for details). 
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At low flow rates (e.g., Vs < lOm/s), the IHS relationship exhibits a linear relation 
r = Xvs, where the friction coefficient A is often used to represent the strength of friction. 
Intrinsically, the friction coefficient is the physically relevant property to characterize the 



interfacial dynamics 



26| . The open squares in Fig. |2t^b) depict the friction coefficients A 
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FIG. 2: (color online) (a) Friction shear stress r versus slip velocity Vg at the water/graphene 
interface with graphene deformation strain e from —10% to 10% (from bottom to top). The solid 
lines represents the fitted relations t/tq = asinh(f5/uo), with the fitting coefficients tq and vq 
summarized in Table [H (b) The friction coefficient A and the slip length Ig as a function of the 
applied strain e. The open square symbols represent the A values calculated from the slopes of 
curves at = 0.0 in (a), and the solid circles are calculated by using the GK relation (Eq. ([2])) in 
our equilibrium MD simulations. The slip length 1^ is calculated by using Ig = ^jl/X [261], where // 
is the viscosity of water. 



(slopes of the t - Vg curves in Fig. [2](a)) as a function of the applied strain e. The key result 
of Fig. [2] is the dramatic effect of the strain on the friction shear stress and the friction 
coefficient. Stretching the graphene layer leads to the increase of the friction stress with a 
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two fold increase of the friction coefficient at the 10% strain. Compressing the graphene 
layer results in a significant reduction in the friction coefficient, e.g., about three fold at 
the —10% strain. In the strain- free state, our calculated friction coefficient is close to the 
value from Falk et al., 1.48 vs. 1.20 (lO^Ns/m^) [27]. The difference can be attributed to 
the different water models and the different LJ potential parameters used to describe the 
carbon-oxygen interactions (possibly yielding different contact angles). Compared to our 
previous MD simulations of water transport inside a double-walled CNT with a diameter 
of 2 nm [qI, the friction coefficient of the graphene channel is ^hig her, i.e., 1.48 vs. 0.3 
(lO^Ns/m^), which is consistent with the conclusion of Falk et al [27|] that A depends on the 
degree of curvature. 

Since the slip length is a widely used quantity to describe the slip flow, we converted our 
friction coefficients A to the slip lengths Is by using /o_= /x/A [26|, where fi is the viscosity of 
water (0.82 mPas for SPC/E water model at 300K 28|]). The results are plot in Fig. [2]^b). 
At the strain-free state (i.e., e = 0), the slip length Is is about 54 nm. Such a large slip 
length in comparison to the channel height of 2 nm implies the plug-like speed profile shown 
in Fig. [D^b). We note that although our calculated slip length is significantly higher than 
the value obtained by Thomas et al. ~ 30 nm [29], the qualitatively form of the agreement is 
reasonable. We believe that the discrepancy arises from the different water model (TIP5P) 
adopted. Overall, however, the remarkable effect shown in these results is the degree to 
which strain engineering of the graphene planes can significantly change the slip length from 
26 nm up to 173 nm. 

The sensitive dependence of the friction coefficient A and the slip length Is on the strain 
(Fig. [2]) suggests that the strain engineering indeed can serve as an effective method of 
controlling the water transport inside the graphene nano-channels. It may also be an im- 
portant factor in our better understanding why there is such scatter in experimental and 
numerical results for slip lengths of water transport inside the CNTs; this can range from 



several nanometers up to several microns lOHlSl, l30|, l3l| . To gain a deep physical insight 



into how the molecular mechanisms affect the drastic changes arising from the imposed 
strain, we investigated the phenomenon further, the results of that research is provided in 
the following sections. 
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B. Molecular Structures of Water inside a Graphene Nano-channel 



It is known that the interfacial molecular structures determine the slip flow over a hy- 
drophobic surface. In this section, we study the structural details of our water/graphene- 
channel system. 

Density profiles of the water molecules across the height of the graphene channels are 
shown in Figure [3t^a). In the strain- free graphene channel, our MD results found almost 
no difference between the densities of the water molecules at rest and at a flow rate of 100 
m/s. The sharp peak indicates the position of the first liquid layer, i.e., 3.25A from the solid 
walls. When different strains (—10% and 10%) are imposed on the graphene planes there is 
no change in the locations of the first density peaks, whereas the peak heights do decrease 
with the strain. The inset of Fig. |3t^a) summarizes the density change with respect to the 
strains. We believe that the reason for the density drop is the weakened attraction between 
the solid walls and the water molecules, arising from the low surface density of the carbon 
atoms caused by the stretching the graphene layer. 

The two-dimension radial density function (RDF) of the first liquid layer was also cal- 
culated. By normalizing the RDF with the mean density, we found that the RDFs almost 
overlap with each other for the differently strained graphene channels, as indicated in Fig. 
int^b) at strains of —10%, 0%, and 10%. This implies that strain engineering on the solid 
walls has a negligible effect on the structures of the first water layer. In comparison, the 
RDF of the bulk water is shown as the dashed line in Fig. [3]^b). The fact that the first peak 
is in almost the position in this case suggests that the average oxygen-oxygen inter-atomic 
distance of the neighbored water molecules are unaffected in the graphene nano-channels, 
i.e.. Too = 0.275 nm. We can thus conclude that the size confinement actually leads to more 
compact water molecules in the first liquid layer than in the bulk water without changing 
the inter-molecular distance. 



The 2-dimensiona^ 



following expression |27|, |32| 



structure factor of the first water layer was calculated by using the 



^ N N ^ / N \ / ^ \ 

j=i 1=1 \j=i J \j=i J 



where Xj represents the position of the jth oxygen atom, is the number of oxygen atoms 
in the first water layer, and q is the 2-dimensional reciprocal lattice vector. The bracket () 
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FIG. 3: (color online) (a) Density profiles of water across the height {z direction) of the graphene 
nano-channel. The sharp density peak next to the solid wall represents the first water layer. For 
graphene layers at zero strain, the density profiles at flow velocity Vg = m/s and Vs = lOOm/s 
have a negligible difference. The density profiles in the stained graphene channels with e = —10% 
(black), 0% (red), 10% (blue) indicate the almost same first peak position but different magnitude. 
Inset shows the density of the first peak as a function of the imposed strains e on graphene. (b) 
Two-dimension radial density function (RDF, normalized by the mean density) of the first water 
layer in three differently strained graphene channels e = —10% (black), 0% (red), 10% (blue). In 
comparison, the RDF of the bulk water is also shown, (c) Two-dimensional structure factor Si{q) 
of the first water layer in the strain free graphene channel exhibits a clear direction independence, 
(d) The radial average of the structure factor in the strained graphene channels with e = —10% 
(black), 0% (red), and 10% (blue). The solid arrows indicate the positions of the reciprocal lattice 
vector \q±\ of the strained graphene planes. 
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denotes an equilibrium ensemble average. Figure |3]^c) shows the structure factor of water 
inside the zero-strained graphene channel. It is clear that the structure factor is almost 
independent of the direction of the q vector. Similar isotropic structure factors are also 
observed in the strained graphene channels. We then averaged the 5'i(q) along all the q 
directions and plot the Si{q) in Fig. El^d) with g = |q| representing the length of the recip- 
rocal lattice vector. The very small differences in the plots of Si{q) confirm our conclusion 
that the strain engineering on graphene walls does not affect the structures of the first water 
layer. However, the strain engineering does change the interfacial commensurateness due to 
the change of the graphene. The solid arrows in Fig. E^d) represent the positions of the 
reciprocal lattice vectors of the differently strained graphene planes |q±|. Here the lattice 
vectors of graphene are connecting the centers of two six-rings in neighbors. With the strains 
from —10% to 10%, it is evident that Sidq^l) increases, which suggests a better degree of 
commensurateness. 

The inset of Fig. H] shows the interfacial potential energy profile over a graphene (e = 0) 
at the position of the first water layer, i.e., zi = 3.25A indicated in Fig. |3]^a). The energy 
profile has a six- fold symmetry. In our analysis, the energy barrier AE{zi) is defined as 
the difference of the potential energy between the maximum points (located on top of the 
carbon atoms) and the minimum points (located in the centre of the six-rings). Figure H] 
depicts a linear relation of the calculated energy barrier AE{zi) with respect to the exerted 
strain. The magnitude of the energy barrier doubles when the strain e changes from —10% 
to 10%. 



C. A Microscopic Understanding on Strain Dependent Friction 

As a dissipation coefficient, the friction coefficient A can be expressed via the Green-Kubo 
(GK) relationship, which relates A to the autocorrelation function of fluctuating pair-wise 



forces at equilibrium 



33| 



A = I {mm)dt (2) 







where F{t) is the total forces in the flow direction exerted on carbon atoms due to the 
interactions with water molecules in our equilibrium MD simulations, fc^ is Boltzmann 
constant, A is the surface area, and T is water temperature. We directly calculated the force 
autocorrelation function (F(t)F(O)) in our equilibrium MD simulations with simulation time 
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FIG. 4: (color online) Interfacial potential energy barriers as a function of the strain applied on the 
graphene walls. The potential energy profile at e = is shown in inset with a six-fold symmetry, 
which is consistent to the atomic structure of grapahene. 

of 5 ns. Here, we should point out that a well-documented difficulty of estimating the GK 



relationships via the equilibrium MD is the finite size of the simulated system, w! 



lich often 



leads to a vanishing of the friction coefficient after a very long time simulation 33|, |3J]. The 
integration should thus have a cutoff time to- A widely adopted method of resolving this is to 
use the onset of a plateau of the integrations as the cutoff to- In our simulations, however, 
it was difficult to locate the plateau in some cases. So we followed the suggestions from 
Ref. 



34 



35[ and chosen to as the first zero of the force autocorrelation function, typically 



this was in the range [Ips, lOps]. The friction coefficients calculated from our equilibrium 
MD simulations (Eq. ([2])) are shown as the solid circles in Figure |2](b). A good agreement 
can be observed in comparison with the NEMD results, which suggests the GK relation is 
applicable in our systems. 
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To gain further insight of the molecular mechanisms, we consider the GK expression (Eq. 
([2])) in more details. The GK relation can be re-expressed as ^3] 

In our equilibrium MD simulations, we find that the de-correlation time Tp = 
/(F(t)F(0))/(F^), weakly depends on the exerted strains on the graphenes. Typically it 
is 100 -160 fs for —10% < e < 10%, which is consistent with the MD results of water flow 



in CNTs (271]. Since the variance of Tp (60%) is one order of magnitude smaller than that 
of A (600%), the main contribution to the variance of A must be from the static root mean 
square (RMS) force (F^). In other words, the change of the friction coefficients A with re- 
spect to the strains should be directly correlated with (F^). Indeed, as shown in Fig. [5](a), 
the good linear relation between A and {F'^)/A suggests that the de- correlation time Tp can 
be approximated as a constant. 

Following Ref. |27|, the RMS force (F^) can be estimated analytically. If we assume 



the main contribution of the friction shear stress arises from the first liquid layer, one can 
approximate the total RMS force as: 



F^) 1 



-pi(5i(q+) + 5i(q_))(goAE)' (4) 



A 2' 

where pi is the density of the first water layer. Si is the 2-dimensinonal structure factor of 
the first water layer, and AE = AE{zi) is the energy barrier. The reciprocal lattice vector 
of the graphene wall is = qo{l / y/S; ±1) with go = 27r/(-\/3acc) and acc the carbon- 
carbon bond length. Figure |5]^b) shows the comparison of the friction coefficients directly 
calculated from our MD simulations (Fig. Mjo)) and those computed by using the analytical 
expressions in Eq. (jl]). The good linear relation clearly validates the form of Eq. (jl]). This is 
an important result because we can quantitatively analyze the effects of density pi, structure 
factor 5*1 (q^), and energy barrier AE on the friction coefficients. 

To quantify the contributions of the energy barrier and the interfacial structures to the 
friction coefficient A, we normalize A by the friction coefficient at zero strain Aq as A/Aq = 
[{qoAE)'^ /{qQAE)l][piSi/{piSi)o]. Figure [6] shows how each of these two terms varies as 
functions of the strain e applied on graphene walls. When compared to Fig. EJ^b), we can 
conclude that: first, the contribution to the change of A mainly comes from the change of 
the energy barrier (Fig. [6](a)); second, in the stretched graphene nano-channels, the increase 
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FIG. 5: (a) Functional dependence of the friction coefficient A versus the static root mean square 
force (F^) (Eq. (b) Functional dependence of the friction coefficient A on the computed 

pi{qQ/S.E)'^ Si (Eq. (jH)). Dashed straight lines are guides to the eye. The good linear relations 
validate the forms of Eqns. ([3]) and 

of commensurateness 5*1 (Fig. Et^d)) and the reduction of the first liquid layer density pi 
(Fig. El^a) inset) cancel out each other, resulting in a small overall contribution (Fig. El^b)); 
third, when the graphene channel under in-plane compression, the decrease of the structural 
factor overwhelms the increase of the density (Fig. EJ^b)), meaning that, quantitatively, the 
contribution from piS*! is about 30% of that from AE towards the effect on the friction 
coefficient (Eq. (HD). 



IV. CONCLUSIONS 



To summarize, our MD simulations show that the friction coefficient A (and the slip 
length Is) of the water transport in the graphenes nanochannels exhibits a highly sensitive 
dependence on the strains imposed on the graphene. The friction coefficients change by 
a factor of 6 when the strain on the graphene wall changes from —10% to 10%. This 
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FIG. 6: (a) The normalized energy barrier of the potential profiles from the graphene solid walls 
as a function of the strains applied, (b) The normalized structural changes as the function of the 
applied strains on the solid walls. 

corresponds to a change in the slip length Ig varies from 26 nm to 173 nm. Our results suggest 
that strain engineering could serve as an effective route to control the water transport inside 
graphene nanochannels. It may also be an important factor in understanding the scatter in 
the reported slip lengths of the water flow inside CNTs in experiments and simulations. The 
molecular mechanisms of the strain effect on the slip flow are also studied. We find that the 
strains on the graphene have relatively small influences on the molecular structure of the first 
water layer, other than on the reduction of the density. Using a simplified analytical model 
based on the Green-Kubo relation [26], we find that the change of energy barrier makes 
the most important contribution to the change of the friction coefficient A, in comparison 
to the effect of the water density and the structural factor (which make about 30% of the 
contribution made by the energy barrier). The quantitative relationship for t — Vs provided 
herein when combined with the physical insights provided on the molecular mechanism will 
be valuable to designers of graphene nano-channels for application in nano-fiuidic devices. 
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